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I.  INTRODUCTION 


Some  of  the  ongoing  concerns  of  the  Weapon  Dynamics  and  Accuracy  Branch  (WDAB)  of 
the  Interior  Ballistics  Division  (IBD)  of  the  Ballistic  Research  Laboratory  (BRL)  are 

(1)  To  realistically  test  instrumentation  and  structures  for  survivability  and  reliability  in 
gun  launch  environments, 

(2)  To  assess  the  damage  to  the  human  ear  due  to  nonclassical  pressure  pulses  such  as 
those  generated  by  multiple  sources  in  confined  volumes, 

(3)  To  establish  economic  screening  and  quality  assurance  procedures  for  shock  testing 
production  items  and 

(4)  To  assess  the  shock  survivability  of  proposed  structural  designs. 


Two  approaches  that  can  be  made  to  these  problems  are 

(1)  To  empirically  duplicate  exactly  the  shock  environment  by  duplicating  the  pulse 
shape,  force  amplitude  and  time  duration  and 

S  *  •  .  •>  w.  4  4*  » 

(2)  To  empirically  duplicate  the  severity  of  the  environment  by  dissimilar  pulses  whose 
severity  has  been  analytically  determined  to  be  relatively  the  same  in  the  frequency 
response  regime  of  the  test  item. 


The  latter  method  has  been  used  successfully  by  scientists  at  BRL  since  the  late  sixties  for 
various  projects.  This  technique  is  dependent  on  the  ability  to  determine  the  shock  spectra  from 
Duhamel’s  integral  for  the  test  and  environmental  pulses,  thus  forming  a  basis  for  comparison  of 
severity.  The  main  advantage  of  the  approach  is  that  it  permits  greater  latitude  in  test  facilities  and 
provides  a  means  of  comparing  different  testing  systems.  It  also  provides  the  designer  with 
information  concerning  the  severity  of  the  environment  and  its  effect  upon  the  structure’s 
components. 

Justification  for  the  method  is  proven  by  consideration  of  several  aspects  of  the  behavior  of 
structures.  First,  all  structural  systems  are  complex  oscillators  and  they  have  responses  in  specific 
frequency  regimes  unique  to  the  particular  systems.  Therefore,  tests  to  determine  survivability  need 
only  reflect  the  energy  and  momentum  existing  in  the  frequency  regimes  of  interest. 

Secondly,  the  relative  response  of  the  structure  is  deterministic  and  repeatable  for  any  given 
pulse  shape. 

Thirdly,  the  response  of  the  structure  can  be  described  by  two  spectra:  the  first  being  the 
primary  spectrum  showing  the  forced  vibration  response  during  the  time  duration  of  the  pulse;  the 
second  being  the  residual  spectrum  showing  the  free  vibration  response  after  the  pulse. 
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n.  BASIC  EQUATION  FOR  DUHAMEL’S  INTEGRAL 


The  spectra  are  obtained  by  calculating  the  response  of  systems  with  discrete  frequencies 
using  the  Duhamel  integral  over  a  time  period  three  times  the  pulse  duration.  The  integral  in  the 
form  used  is  shown  in  Eq.  (1)  . 


Ar  =  I  x0w* - —  J  A  (tv)  sin(w„fv)  dtv  I  cos(w„j) 


+  I  X(/Jin  +  —  J  A(ty)  cos(w„ tv )  dty  I  sin(w„  j  ) 


(1) 


where 


Ar  =  response  acceleration 
x0  =  initial  displacement  of  the  system 
x0  =  initial  velocity  of  the  system 
ojn  =  natural  frequency  of  the  system 
A  (tv)  =  instantaneousaccelerationof  the  pulse  at  time  tv 
tv  =  time  of  observation 
s  =  period  of  the  system  of  frequency  u„. 


A  numerical  approach  was  necessary  because  the  parameter  that  describes  the  motion 
(A  ( ))  is  generally  an  unknown  function  that  is  approximated  by  a  set  of  ordered  pairs  taken  from 
a  digitized  set  of  experimental  firing  data. 


If  we  assume  that  the  pulse  starts  with  zero  displacement  and  velocity,  then  Eq.  (1)  can  be 
simplified  to  that  shown  in  Eq.  (2)  . 


Ar  = 


sin  (w,i) 


/  A(ty)  COS(Wn  ty)  dty  ~ 


cos(< jns) 


/  A  (ty)  S\n(U)n  ty)  dty 


(2) 


Since  the  sum  of  the  integral  equals  the  integral  of  the  sum,  Eq.  (2)  can  be  expanded  so  that 
we  do  not  integrate  from  zero  for  each  point  of  the  time  history .  The  natural  frequency  un  of  the 
system  in  Eq.  (3)  is  replaced  by  the  damped  frequency  wf.  This  damped  frequency  is  calculated  by 
multiplying  the  natural  frequency  by  \/l-62  ,  where  6  is  the  damping  factor.  The  expanded  Eq. 
(2)  with  the  damped  frequency  is  shown  in  Eq.  (3)  . 
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(3) 


Ar  =  —  I  sin  (w*  )  \  f  A  (ty)  cos (ws  tv)  dtv  +  J  A  (fv )  cos(w«  tv )  dtv 

(j)^  L.  I  *.tn  ■  Cj 


<-H  \*s  0 


cos(o;4  )  ]  I  A  (ty)  sm(w6  tv)  dtv  4-  /  A  (tv)  sm(u6  tv)  dtv 

i-H 


] 


Because  it  is  not  necessary  to  integrate  from  zero  for  each  point  of  time  history,  Eq.  (3) 
provides  considerable  savings  in  central  processing  unit  (cpu)  time.  Hence,  Eq.  (3)  was  used  in  the 
program  instead  of  Eq.  (2).  This  equation  produces  the  primary  system  response  for  a  given 
frequency. 

A  damping  coefficient,  which  takes  into  account  the  viscous  loss  due  to  velocity  in  the 
system,  has  been  incorporated  into  Eq.  (3).  This  coefficient  was  determined  from  Eq.  (4). 

Dampp  =  (4) 

where  Dampp  =  damping  coefficient  for  primary  system  response 
8  =  damping  factor 
wn  =  natural  frequency  of  the  system 
s  =  period  of  the  system  of  frequency  w„ . 

The  residual  system  response  is  obtained  analytically  from  Eq.  (5). 

R,  =  Arz  cos  -z  ))  +  A  'rz  sin(u >^t  -z  ))  (5) 

where  R,  =  residual  response  acceleration  at  time  z 
Arz  =  response  acceleration  at  time  z 
A  'rz  =  response  velocity  at  timez 
t  =  time  after  pulse  ends 
ws  =  damped  frequency  of  the  system 
z  =  pulse  duration. 
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As  in  the  primary  system  response,  the  natural  frequency  in  Eq.  (5)  has  been  replaced  by  a 
damped  frequency.  Eq.  (6)  shows  the  resulting  formulation. 


A  V 

R  't  =Arz  e ("**'*'■*))  cos (w^t-z))  H - - — — - sin(w^(f  -z )) 


(6) 


(5). 


Where  the  definitions  of  all  variables  are  the  same  as  those  previously  listed  in  Eqs.  (4)  and 


III.  THE  COMPUTER  PROGRAM 


The  BRL  mainframe  computer  during  the  past  few  years  has  been  a  Control  Data 
Corporation  (CDC)  CYBER  7600(MFZ).  The  two  front-end  machines  were  a  CYBER  750(MFA) 
and  an  825(MFB).  Two  new  BRL  supercomputers,  the  CRAY  XMP/48  and  the  CRAY  2  have 
been  installed  at  the  BRL  site.  There  are  various  minicomputers  and  microcomputers  within  the 
divisions  of  BRL.  The  IBD  has  several  Hewlett-Packard  (HP)  microcomputers,  namely  the 
HP9845C,  the  HP9836C  and  the  HP9000.  Also,  within  the  IBD  is  an  HP1000F  minicomputer. 

A  computer  program  was  written  to  numerically  evaluate  Duhamel’s  integral.  The  numerical 
integration  scheme  used  in  this  program  is  Simpson’s  integration.  The  reciprocal  of  a  given 
frequency  is  used  as  the  period  of  the  function.  Additionally,  the  step  size  for  Simpson’s  integration 
is  less  than  one-eighth  the  period  to  attain  reasonable  accuracy  using  this  numerical  integration 
scheme.  This  step  size  scheme  was  selected  so  that  the  program  would  use  fewer  calculations  in 
the  lower  frequencies. 

The  natural  frequency  w  of  the  system,  given  by  u  =  2%f  where  /  is  the  inputted  frequency 
(in  Hertz),  is  an  integral  part  of  the  function  being  integrated  as  illustrated  in  Eqs.  (1)  through  (6). 

The  Duhamel  integral  is  currently  coded  in  both  FORTRAN  and  BASIC  programming 
languages  and  has  been  run  on  the  Hewlett-Packard  (HP)  984$C,  9836C,  CYBER  750,  825  and 
7600  and  the  CRAY  XMP/48.  The  code  can  also  be  used  on  the  HP1000,  HP9000  and  other 
machines  which  support  FORTRAN  and/or  BASIC.  The  FORTRAN  version  is  currently  being 
prepared,  along  with  its  job  control  language,  to  be  run  on  the  CRAY  2  machine. 

The  program  is  fully  documented  and  is  generic  for  both  damped  and  undamped  systems. 
The  FORTRAN  code  along  with  its  associated  job  control  language’ for  CRAY  XMP/48  is  listed  in 
Appendix  A. 

The  functions  and  interrelationships  of  files  run  on  CRAY  XMP/48  are  found  in  Figures  1 
and  2.  Figure  3  is  a  general  flow  diagram  of  the  basic  program.  ^ 
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•  Job  Control  Longuoge 

Figure  1.  CRAY  XMP/48  Integration  Program  Sequence 


Job  Control  Language 


Figure  2.  CRAY  XMP/48  3-D  Plot  Program  Sequence 


11 


(  start") 


Read  or  Generote  Pulse  Doto 
Acceleration  vs  Time 


<D-» 


Initiolize  Porameters 
ond 

Read  Frequencies 


Interpolate  for 
Equally-spoced  Points 


P  Integrate  Duhamel’s  Equotion 
far  Points  on  Time-History  Curve 


N  / 


Compute  Moximum  Time  History 
for  Primary  ond  Residuol  Ports 


Frequencii 
_  Seen  Read 


©  lncf<owl  t/RS 


Time 


[  Yes 

Plot  Spectra 

* 

£ _ „ 

CstoD 


Figure  3.  General  Flow  Diagram  of  the  Basic  Integration  Program 


A.  Input  to  the  Computer  Program 


The  first  step  in  establishing  test  criteria  is  to  obtain  response  spectra  of  the  shock  pulse 
representing  the  environment.  Thus  far,  the  half-sine  pulse,  test  machine  pulses,  interior  ballistic 
pulses,  idealized  free-space  blast,  modified  free-space  blast  and  measured  confined  blast  pulses  have 
been  used  as  representations  of  the  environment.  Figures  4  through  8  are  examples  of  these 
different  types  of  shock  pulses. 

These  types  of  shock  pulses  are  used  as  input  to  the  program.  Figure  4,  a  half-sine  pulse,  is 
an  example  of  a  code-generated  curve  which  was  generated  in  the  Duhamel  code.  The  idealized 
free-space  blast  and  modified  blast  pulses,  Figures  5  and  6,  are  curves  that  were  generated  in  an 
auxiliary  program  for  utilization  in  the  Duhamel  code.  Figures  7  and  8  are  examples  of  actual 
experimental  data  curves  taken  from  firing  data. 

The  amplitude  of  the  input  data  curves  is  normalized  to  one,  thus  eliminating  the  need  for 
dimensions.  The  frequency  axis  is  normalized  by  multiplying  the  frequency  by  the  pulse  duration. 
This  yields  cycles.  Depending  on  the  program  option  selected,  the  frequency  characteristicscan  be 
inputted  from  the  keyboard,  read  from  the  data  statement  within  the  program  or  included  as  a  part 
of  the  input  file.  The  frequency  range  is  generally  between  .1  and  10000  cycles.  The  damping 
characteristics  are  inputted  from  the  keyboard  or  from  a  data  file  and  range  from  zero  to  one. 
Figure  9  is  an  example  of  input  data.  Table  1  defines  the  input  variables  and  lists  the  formats. 
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AMPLITUDE 


Figure  4.  Program  Input  -  Half-Sine  Pulse 


Figure  5.  Program  Input  -  Idealized  Free-Space  Blast  Pulse 
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PRESSURE.  CTS  RMPLITUDE 


TIME 

Figure  6.  Program  Input  -  Modified  Blast  Pulse  (Positive  Portion) 


si 


TIME,  MS 

Figure  7.  Program  Input  -  Pressure  Pulse  Tera  Socorro  ID  24 
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Figure  8.  Program  Input  -  Acceleration  Pulse  105-mm  Gun  M68  M456A2 


Card 

Number 

■  '  Data 

1 

201  0  3.14159  1.0  10000.  100.  0  1  1  1  0 

2 

105 MM  GUN,  M68  M456A2 

3 

ACC-TIME  HISTORY  OF  M456A2 

4 

TIME,  MSEC  ACCELERATION,  G’S 

7a 

•000 

10 

SHOCK  SPECTRA  M456A2 

11 

FRQ  *  PULSE  DURATION  STATIC  ACCELERATION 

12 

.1  10000.0  1000.  0.0  5.0  1. 

Figure  9.  Sample  FORTRAN  Card  Image  Input  For  Duhamel  Integral  Code 
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TABLE  1.  DEFINITION  AND  FORMAT  OF  INPUT  VARIABLES  TO 
DUHAMEL  INTEGRAL  CODE 


Record 

Number 

Variable 

Name 

Column 

Numbers 

Format 

Description 

1 

NOC 

1-5 

15 

Number  of  input  points. 

IREAD 

6-10 

15 

Input  option,  if  IREAD  =  1,  then 
read  the  input  data  curves  from 
TAPE  1,  otherwise,  generate  curve 
in  program. 

TCAL 

11-20 

F10 

Time  duration  of  pulse. 

ACAL 

21-30 

F10 

Maximum  amplitude  of  input  pulse. 

FMAX 

31-40 

F10 

Maximum  frequency  times  pulse  duration. 

TS 

41-45 

F5 

Number  of  points  to  be  generated 
for  each  time  history. 

ILOG 

46-50 

15 

Logarithmic  option,  if  ILOG  =  1, 
then  calculate  and  plot  logarithmic 
values  of  frequency. 

IRES 

51-55 

15 

Residual  option,  if  IRES  =  1,  then 
plot  only  residual  spectrum. 

IFREQ 

56-60 

15 

Frequency  option,  if  IFREQ  =  1, 
then  convert  frequency  times  pulse 
duration  to  frequency. 

ISPEC 

61-65 

15 

Spectra  option,  if  ISPEC  =  1,  then 
plot  both  primary  and  residual 
spectra. 

I3D 

66-70 

15 

3-D  option,  if  I3D  =  1,  then  plot  a 

3-D  spectra  surface. 

IPLOR 

71-72 

12 

Plot  option,  if  IPLOR  =  1,  then 
plot  the  original  input  pulse, 
unnormalized. 

IPLNO 

73-74 

12 

Plot  option,  if  IPLNO  =  1,  then 
plot  the  original  input  pulse, 
normalized. 
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TABLE  1.  DEFINITION  AND  FORMAT  OF  INPUT  VARIABLES  TO 
DUHAMEL  INTEGRAL  CODE(cont.) 


Record 

Number 

Variable 

Name 

Column 

Numbers 

Format 

Description 

2 

PGTIT 

1-80 

10A8 

Page  title  to  be  printed  on 
tabulated  output. 

3 

i 

TIT4 

1-80 

10A8 

Plot  title  for  unnormalized  input 
pulse. 

4 

LXNAME 

1-24 

3A8 

Label  for  x  axis  for  plot  of 
unnormalized  input  pulse. 

LYNAME 

25-48 

3A8 

Label  for  y  axis  for  plot  of 
unnormalized  input  pulse. 

5 

TIT4 

1-80 

10A8 

Plot  title  for  normalized  pulse. 

6 

LXNAME 

1-24 

3A8 

Label  for  x  axis  for  plot  of 
normalized  input  pulse. 

LYNAME 

25-48 

3A8 

Label  for  y  axis  for  plot  of 
normalized  input  pulse. 

7 

FT 

1-10 

F10 

Frequency  times  pulse  duration. 

DAMP 

11-20 

F10 

Damping  characteristic. 

IPLOT 

21-25 

15 

Plot  option,  if  IPLOT  =  1,  then 
plot  the  time  history  of  the 
individual  frequency. 

7a 

DAMP 

1-10 

F10 

If  a  spectra  is  being  plotted,  then 
the  damping  characteristic  is  the 
only  variable  on  Input  Card  7  that 
needs  to  be  read.  It  only  needs  to 
be  read  once,  not  for  each 
frequency. 

8 

TIT4 

1-80  . 

10A8 

Plot  title  for  plot  of  time  history  of 
an  individual  frequency. 

9 

LXNAME 

1-24 

3A8 

Label  for  x  axis  for  plot  of  time 
history  of  an  individual  frequency. 
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TABLE  1.  DEFINITION  AND  FORMAT  OF  INPUT  VARIABLES  TO 
DUHAMEL  INTEGRAL  CODE(cont.) 


Record 

Number 

Variable 

Name 

Column 

Numbers 

Format 

Description 

9 

LYNAME 

25-48 

3A8 

Label  for  y  axis  for  plot  of  time 
history  of  an  individual  frequency. 

10 

TIT4 

1-80 

10A8 

Plot  title  for  plot  of  residual  and/or 
primary  spectra. 

11 

LXNAME 

1-24 

3A8 

Label  for  x  axis  for  plot  of  residual 
and/or  primary  spectra. 

LYNAME 

25-48 

3A8 

Label  for  y  axis  for  plot  of  residual 
and/or  primary  spectra.  * 

12 

XORIG 

1-10 

E10.3 

I 

Minimum  x-scale  value  for  plot  of 
residual  and/or  primary  spectra. 

XMAX 

11-20 

E10.3 

Maximum  x-scale  value  for  plot  of 
residual  and/or  primary  spectra. 

XSTP 

21-30 

E10.3 

Data  units  per  x-scale  plot  inch  for 
plot  of  residual  and/or  primary 
spectra. 

YORIG 

31-40 

E10.3 

Minimum  y-scale  value  for  plot  of 
residual  and/or  primary  spectra. 

YMAX 

41-50 

E10.3 

Maximum  y-scale  value  for  plot  of 
residual  and/or  primary  spectra. 

YSTP 

51-60 

E10.3 

Data  units  per  y-scale  plot  inch  for 
plot  of  residual  and/or  primary 
spectra. 

B.  Output  From  the  Computer  Program 

Two  spectra  are  obtained  as  output  from  the  Duhamel  code. 

The  primary  spectrum  which  is  the  maximum  positive  or  negative 
relative  acceleration  achieved  during  the  pulse. 

The  residual  spectrum  which  is  the  maximum  positive  or  negative 
relative  acceleration  achieved  after  the  pulse. 


Both  spectra  are  plotted  versus  frequency. 
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The  tabulated  output  of  the  Duhamel  code  includes  the  input  in  its  original  form,  the  input  in 
its  normalized  form  and  the  individual  responses  for  the  selected  frequencies. 

The  plotted  output  is  a  choice  of  any  or  all  of  the  following: 

(a)  the  input  curve  in  original  form, 

(b)  the  input  curve  in  normalized  form, 

(c)  individual  frequency  response, 

(d)  a  spectra  of  the  primary  and  residual  system  responses  to  the  given 

frequencies,  and/or 

(e)  a  three-dimensionalmappingof  the  primary  and  residual  system  responses 

showing  relative  acceleration,  time  and  frequencies. 

C.  Plot  Routines 

On  the  Hewlett-Packard  machines,  there  is  a  separate  plot  routine  for  each  of  the  plot 
options.  The  data  are  stored  on  disc  and  then  read  in  by  the  selected  plot  program.  In  the  case  of 
the  3-D  option,  on  all  of  the  machines,  the  data  must  be  stored  for  each  system  response 
individually.  The  data  are  then  read  by  the  3-D  routine,  geometrically  rotated  and  subsequently 
plotted.  The  rotation  is  done  to  provide  an  optimum  view  of  each  of  the  individual  system 
responses.  The  FORTRAN  code  for  the  3-D  plot  and  the  associated  job  control  language  for  the 
CRAY  XMP/48  are  listed  in  Appendix  B.  Figure  10  is  an  example  of  actual  input  data  while 
Table  2  provides  a  listing,  including  descriptions  and  formats,  of  the  input  variables  for  the  3-D 
spectra  plot. 

On  the  Hewlett-Packard  machines,  the  HP  graphics  routines  are  used,  while  on  the  CRAY 
XMP/48,  the  graphics  are  done  using  the  commercial  plotting  package  DISSPLA.  Version  9  of 
DISSPLA  is  now  resident  on  the  CRAY  XMP/48.  The  plots  created  using  DISSPLA  can  be 
displayed  on  a  terminal  screen,  plotted  on  a  printer  attached  to  an  interactive  terminal  and/or 
transferred  to  the  CALCOMP  plotter  at  the  BRL  central  computer  site  for  hardcopy. 


Card 

Number 

Data 

1 

50  2  "  - - 

2 

105MM  GUN,  M68  M456A2  ~  “ 

3 

0.  2.5  .3125  -4. 

1.5  .9167 

Figure  10.  Sample  FORTRAN  Card  Image  Input  for  Duhamel  3-D  Spectra  Plot 


19 


Table  2.  DEFINITION  AND  FORMAT  OF  INPUT  VARIABLES  FOR 
DUHAMEL  3-D  SPECTRA  PLOT 


Record 

Number 

Variable 

Name 

Column 

Numbers 

Format 

Description 

1 

NFILES 

1-5 

15 

Number  of  response  files  to  be  plotted. 

MAXT 

6-10 

15 

Maximum  time  of  response  plot. 

1  =  Primary  response  only 

2  =  Primary  and  partial 

residual  responses 

3  =  Primary  and  full  residual 

responses 

2 

TIT4 

1-80 

10A8 

Plot  title 

3 

XORIG 

1-10 

El  0.3 

Minimum  x-scale  value 

XMAX 

11-20 

E10.3 

Maximum  x-scale  value 

XSTP 

21-30 

El  0.3 

Data  units  per  x-scale  plot  inch 

YORIG 

31-40 

E10.3 

Minimum  y-scale  value 

YMAX 

41-50 

El  0.3 

Maximum  y-scale  value 

YSTP 

51-60 

El  0.3 

Data  units  per  y-scale  plot  inch 

Figures  11,  12  and  13  show  the  response  of  systems  with  30  cycles  normalized  frequency  and 
.05  damping  to  the  pulses  in  Figures  5  through  7,  respectively.  Figures  11  through  13  are  individual 
system  responses  versus  time  for  a  single  frequency.  Figure  14  is  an  overlay  of  Figures  11,  12  and 
13. 


Figure  15  is  the  response  spectra  of  systems  with  .1  to  10000  cycles  normalized  frequencies 
and  .05  damping  applied  to  the  half-sine  pulse,  Figure  6.  Figures  16  through  20  are  response 
spectra  of  systems  with  .1  to  10000  cycles  normalized  frequencies  and  no  damping  applied  to 
Figures  4  through  8,  respectively.  Figures  11  through  14  show  both  primary  and  residual  response 
spectra  versus  frequency  plotted  on  a  semi-log  plot. 

In  actuality,  these  represent  normalized  response  curves  where  the  abscissa  is  frequency  times 
pulse  duration  (cycles)  and  the  ordinate  is  in  relative  acceleration  (dimensionless). 
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RELATIVE  AMPLITUDE  RELATIVE  AMPLITUDE 


Figure  11.  Response  Curve  (Frequency  =  30,Damping  =  .05)  Idealized 

Free-Space  Blast  Pulse 


Figure  12.  Response  Curve  (Frequency  =  30,  Damping =.05)  Modified  Blast  Pulse 

(Positive  Portion) 
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RELRTIVE  AMPLITUDE 


Figure  13.  Response  Curve  (Frequency  =  30,  Damping  =  .05)  Pressure  Pulse 

Tera  Socorro  ID  24 


Figure  14.  Comparison  Response  Curves  (Frequency  -  30,  Damping  -  .05) 
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NORMALIZED  SHOCK  RESPONSE  SPECTRA 


PRIMARY 
•  RESIDUAL 


Figure  15.  Normalized  Shock  Response  Spectra  -  Half-Sine  Pulse  (Damping  =  .05) 

NORMALIZED  SHOCK  RESPONSE  SPECTRA 


PRIMRRY 
•  RESIDUAL 


Figure  16.  Normalized  Shock  Response  Spectra  -  Half-Sine  Pulse  (Undamped) 
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NORMALIZED  SHOCK  RESPONSE  SPECTRA 


-  PRIMARY 

- « -  RESIDUAL 


Figure  17.  Normalized  Shock  Response  Spectra  -  Idealized  Free-Space  Blast 

Pulse 

NORMRLIZED  SHOCK  RESPONSE  SPECTRA 


PRIMARY 

RESIDUAL 


Figure  18.  Normalized  Shock  Response  Spectra  -  Modified  Blast  Pulse 

(Positive  Portion) 
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RELATIVE  AMPLITUDE  RELATIVE  AMPLITUDE 


NORMALIZED  SHOCK  RESPONSE  SPECTRA 


PRIMARY 

RESIDUAL 


FREQUENCY  TIMES  PULSE  DURATION 


Figure  19.  Normalized  Shock  Response  Spectra  -  Pressure  Pulse 
Tera  Socorro  ID  24 

NORMALIZED  SHOCK  RESPONSE  SPECTRA 


PRIMARY 

RESIDUAL 


Figure  20.  Normalized  Shock  Response  Spectra  -  Acceleration  Pulse 
105-mm  Gun  M68  M456A2 
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Figures  21  through  26  are  3-D  mappings  of  the  response  spectra  shown  in  Figures  15  through 
20,  respectively.  The  three  dimensions  are  relative  acceleration,  time  and  frequency. 


Figure  21.  3-D  Shock  Response-  Half-Sine  Pulse  (Damping  =  .05) 


* 


Figure  22.  3-D  Shock  Response-  Half-Sine  Pulse  (Undamped) 
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Kggg 

i 

Figure 23.  3-D  Shock  Response-  Idealized Free-Space Blast  Pulse 


■! 


Figure  24.  3-D  Shock  Response-  Modified  Blast  Pulse  (Positive  Portion) 
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Figure  25.  3-D  Shock  Response-  Pressure  Pulse  Tera  Socorro  ID  24 


Figure  26.  3-D  Shock  Response-  Acceleration  Pulse  105-mm  Gun  M68  M456A2 
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IV.  SUMMARY 


The  Duhamel’s  integral  computer  code  was  developed  to  analyze  the  transient  response  of 
test  instrumentation  and  structures  and  to  assess  damage  to  the  human  ear  caused  by  nonclassical 
pressure  pulses.  The  code  has  been  used  to  assess  the  relative  severity  of  laboratory  test 
environments  with  respect  to  gun  environments.  The  program  has  been  installed  on  a  mainframe,  a 
front-end  machine  and  several  microcomputers  in  BASIC  and  FORTRAN  programming  languages 
so  as  to  make  it  more  accessible  to  the  user. 

The  Duhamel’s  integral  package  is  now  being  run  on  the  CRAY  2  supercomputer  which, 
along  with  the  CRAY  XMP/48,  has  replaced  the  CYBER  7600  at  the  BRL.  This  code  will  be 
incorporated  into  the  general  structural  and  analysis  scheme  of  the  IBD.  In  addition,  the  program 
will  be  extended  to  compound  response  descriptions  which  include  more  than  one  frequency  and 
other  nonlinear  forms,  such  as  hardening  and  softening  stiffness  systems. 
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APPENDIX  A 

LISTING  OF  FORTRAN  PROGRAM  TO  INTEGRATE  DUHAMEL’S  INTEGRAL 
ALONG  WITH  ASSOCIATED  JOB  CONTROL  LANGUAGE  FOR  CRAY  XMP/48 
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APPENDIX  A 


LISTING  OF  FORTRAN  PROGRAMS  TO  INTEGRATE  DUHAMEL’S  INTEGRAL 
ALONG  WITH  ASSOCIATED  JOB  CONTROL  LANGUAGE  FOR  CRAY  XMP/48 

I.  Job  Control  Language  for  CRAY  XMP/48 

JOB,JN = jobname  ,T=  150,MFL=350Q00.  [  ' 

ACCOUNT yfii.C= accountnumber  >\JS=usemame , UPW =userpassword. 

ACCESS,  DN  =  DISLIB,PDN = DISLJB,ID  =  DISSPLA.OWN = SYSTEM. 

ACCESS, DN  =  D  VSD,PDN  =  D  VSD,ID  =DlSSPLA,0  WN  =  SYSTEM. 

ACCESS, DN  =  INTLIB,PDN  =  INTLIB,ID  =  DISSPLA,OWN  =  SYSTEM. 

FETCH, DN = infile  ,TEXT = ’CHARG Epcctno ,  pn  .GET,  infile  =infile  .CTASK.’. 

(  A  .* '  * '*'  ,7  '  *  ■  ?  i 

.  ■  ;  17  .  '  "  •  '  • 

where  infile  =  name  of  input  file  which  contains  program  options 

and  plot  titles 

pn  =  first  two  letters  of  the  account  number 
ASSIGN, DN=infile  A  =  FT03. 

FETCH,DN=a  c/i/e  ,TEXT  =  ’CHARGE,  acctno ,  pn  .GET  pc  file  =acfile  .CTASK.’. 

where  acfile  =  name  of  permanent  file  residing  on  front  end  which 
contains  data  to  be  used  as  shock  pulse 

ASSIGN, DN = acfile  A  =  FT01 . 

ASSIGN, DN  =  TAPE2A  =  FT02. 

CFT. 

SEGLDR,CMD  =  ’LIB  =  DISLIB,INTLIB’,GO. 

DISPOSE, DN =META,DC= ST,DF  =  BB,~ 

TEXT  =  ’CHARGE  =  acctno ,  pn  .CTASK.R  EPLACE,META  =  pfn  .’. 

where  pfn  =  name  of  output  file  in  which  the  DISSPLA  plot  file 
is  stored  on  the  front  end. 

DISPOSE, DN  =  TAPE2,DC  =  ST,DF = BB,^ 

TEXT  =  ’CHARGE = acctno ,  pn  .CTASK.REPLACE,TAPE2  =  3dfile .’. 

where  3d  file  =  name  of  output  file  that  will  contain  data  to 
be  used  in  calculations  for  3-D  plot. 
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on  non 


II.  Listing  of  FORTRAN  File 


C 

C 

C 

C 

C 

C 

C 

C 


C 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PROGRAM  SPECTRM(INPUT, OUTPUT, TAPE1,TAPE5,TAPE6 = OUTPUT, 
*META,TAPE10 = OUTPUT, TAPE2) 

THIS  PROGRAM  NUMERICALLY  INTEGRATES  DUHAMEL’S  INTEGRAL  TO 
OBTAIN  A  PRIMARY  AND  RESIDUAL  SHOCK  SPECTRA.  THE  INPUT 
CONSISTS  OF  A  SET  OF  ORDERED  PAIRS  THAT  DESCRIBE  THE 
ACCELERATION  OF  A  SYSTEM  AND  A  SET  OF  FREQUENCIES  OVER 
WHICH  THE  SYSTEM  IS  TO  BE  EVALUATED.  THIS  PROGRAM  WILL  ALSO 
PLOT  THE  PRIMARY  AND  RESIDUAL  SHOCK  TIME  HISTORY  OF  A  GIVEN 
FREQUENCY. 


DIMENSION  RMAX(500),FWN(500),FR(500),MSAVE(75),XPLT(7000) 
DIMENSION  PFT(200),PXP(200),PXR(200),PGTIT(10),YPLT(7000) 
DIMENSION  XA(3500),XB(3500),TA(3500),TB(3500),TIT4(4) 
DIMENSION  U(210),A(210),V(1040),VP(1040),DUM(1),FTR(61) 

REAL  NX,  NXP,  MAXP,  MAXR 

DATA  PI2, PI/6.283185307, 3.14159265358/ 

FREQUENCY  INPUT  DATA 


DATA  FTR/.l,. 133, .167, .2, .25, .3, .35, .4, .5, .6, .7, .8, .9,1., 1.333, 


*1.667, 2., 2.5, 3., 3.5, 4., 5., 6., 7., 8., 9., 10., 13.333, 16.667, 20., 25., 

*30., 35., 40., 50., 60., 70., 80., 90., 100., 133.333, 166.667, 200., 250., 
*300., 400., 500., 600., 700., 800,900., 1000., 2000., 3000., 4000,5000., 
*6000,7000,8000,9000,10000./ 


CALL  COMPRS 
CALL  SETDEV(10,6) 

ACCELERATION  AND  TIME  ARE  NORMALIZED. 

NOC  =  NO.  OF  INPUT  POINTS. 

IF  IREAD  =  1,  THEN  READ  INPUT  DATA  FROM  TAPE1 

OTHERWISE,  GENERATE  DATA  IN  PROGRAM. 

TCAL  =  TIME  DURATION  OF  PULSE 
ACAL  =  MAXIMUM  ACCELERATION  OF  PULSE 
FMAX  =  MAXIMUM  FREQUENCY  TIMES  PULSE  DURATION 
TS  =  THE  NUMBER  OF  POINTS  GENERATED  FOR  EACH  TIME  HISTORY. 
IF  I  LOG  =  1,  CALCULATE  AND  PLOT  LOG  VALUES  OF  FREQUENCY 
IF  IRES  =  1,  PLOT  ONLY  RESIDUAL  SPECTRUM 
IF  IFREQ  =  1,  CONVERT  FREQUENCY  TIMES  PULSE 
DURATION  TO  FREQUENCY. 

IF  ISPEC  =  1,  PLOT  PRIMARY  AND  RESIDUAL  SPECTRA 
IF  I3D  =  1,  PLOT  A  3-D  SPECTRUM  SURFACE. 

IF  IPLOR  =  1,  PLOT  THE  ORIGINAL  PULSE,  UNNORMALIZED. 

IF  IPLNO  =  1,  PLOT  THE  ORIGINAL  PULSE,  NORMALIZED. 

34 


SPEC  1 
SPEC  2 
SPEC  3 
SPEC  4 
SPEC  5 
SPEC  6 
SPEC  7 
SPEC  8 
SPEC  9 
SPEC  10 
SPEC  11 
SPEC  12 
SPEC  13 
SPEC  14 
SPEC  15 
SPEC  16 
SPEC  17 
SPEC  18 
SPEC  19 
SPEC  20 
SPEC  21 
SPEC  22 
SPEC  23 
SPEC  24 
SPEC  25 
SPEC  26 
SPEC  27 
SPEC  28 
SPEC  29 
SPEC  30 
SPEC  31 
SPEC  32 
SPEC  33 
SPEC  34 
SPEC  35 
SPEC  36 
SPEC  37 
SPEC  38 
SPEC  39 
SPEC  40 
SPEC  41 
SPEC  42 
SPEC  43 
SPEC  44 
SPEC  45 
SPEC  46 
SPEC  47 


c 

SPEC  48 

c 

SET  UP  PLOT  VARIABLES 

SPEC  49 

c 

SPEC  50 

ISC=0 

SPEC  51 

ITYPE=0 

SPEC  52 

N4=0 

SPEC  53 

ICOLOR=0 

SPEC  54 

NCURV = 1 

SPEC  55 

c 

SPEC  56 

c 

READ  AND  PRINT  PROGRAM  CONTROL  VARIABLES 

SPEC  57 

c 

SPEC  58 

10 

READ(5,1000)  NOC,IREAD,TCALrACAL,FMAX,TS,ILOG,IRES, 

SPEC  59 

*IFREQ,ISPEC,I3D,IPLOR,IPLNO 

SPEC  60 

IF(EOF(5).EQ.O.)  GO  TO  20 

SPEC  61 

CALL  DONEPL 

SPEC  62 

STOP 

SPEC  63 

20 

WRITE(6,1000)  NOC,IREAD,TCAL,ACAL,FMAX)TS,ILOG,IRES, 

SPEC  64 

*IFREQ,ISPEC,I3D,IPLOR,IPLNO 

SPEC  65 

c 

SPEC  66 

c 

READ  IN  TITLE  AND  SET  UP  OTHER  PLOT  VARIABLES 

SPEC  67 

c 

SPEC  68 

READ(5,1100)  PGTIT 

SPEC  69 

WRITE(6,1200)  PGTIT 

SPEC  70 

IF(IFREQ.EQ.l)  FMAX  =  FMAX /TCAL 

SPEC  71  . 

c 

SPEC  72 

c 

DO-LOOP  TO  READ  AND  NORMALIZE  ORDERED  PAIRS  (U,A) 

SPEC  73 

c 

SPEC  74 

IF(IREAD.NE.l)  GO  TO  30 

SPEC  75 

READ(1,1600)  (U(I)^(I),I  =  l,NOC) 

SPEC  76 

GO  TO  50 

SPEC  77 

30 

TIME=0. 

SPEC  78 

DO  40  I  =  l,NOC 

SPEC  79 

ACC = SIN(TIME) 

SPEC  80 

U(I)=TIME 

SPEC  81 

A(I)=ACC 

SPEC  82 

TIME = TIME  +  PI /(NOC-1) 

SPEC  83 

40 

CONTINUE 

SPEC  84 

50 

WRITE(6,1300) 

SPEC  85 

WRITE(6,1400) 

SPEC  86 

WRITE(6,1500)  (U(I)rA(I),I  =  l,NOC) 

SPEC  87 

C 

SPEC  88 

C 

PLOT  ORIGINAL  INPUT  PULSE 

SPEC  89 

C 

SPEC  90 

IF(IPLOR.NE.l)  GO  TO  60 

SPEC  91 

READ(5,1100)  TIT4 

SPEC  92 

LST=1 

SPEC  93 

LSP  =  NOC 

SPEC  94 

CALL  DATA4(NCURV,UrA,DUM,NOC, ISC, ITYPE,LST,LSP,TIT4, 

SPEC  95 

*ICOLOR) 

SPEC  96 

TIME=0. 

SPEC  97 

35 


non 


c 

c 

c 

60 


70 

c 

c 

c 


c 

c 

c 


80 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NORMALIZE  INPUT  PULSE 

DO  70  I=l,NOC 
U(I) = U  (I)  /TCAL 
A  (I)  =  A(I)/ACAL 
UMAX=U(I) 

CONTINUE 

PRINT  NORMALIZED  ACCELERATION  PULSE 

WRITE(6,1700) 

WRITE(6,1300) 

WRITE(6,1500)  (U(I)rA(I),I  =  l,NOC) 

PLOT  PULSE  CURVE 

IF(IPLNO.NE.l)  GO  TO  80 

READ(5,1100)  TIT4 

LST=1 

LSP=NOC 

NCURV = 1 

CALL  DATA4(NCURV,U^,DUM,NOC, ISC, ITYPE,LST,LSP,TIT4 
♦ICOLOR) 

TIME  AND  ACCELERATION  MUST  START  AT  0. 


A(1)=0. 

U(1)=0. 

S  =  RANGE  OVER  WHICH  DUHAMEL’S  INTEGRAL  IS  INTEGRATED 
TO  OBTAIN  A  POINT  OF  TIME  HISTORY  FOR  PRIMARY  SHOCK 
SPECTRUM.  PRIMARY  SHOCK  SPECTRUM  IS  GENERATED  AS  XA  VS  S. 


s=o 


T  IS  USED  TO  GENERATE  EQUALLY  SPACED  POINTS  FOR 
SIMPSON’S  INTEGRATION 


T  =  0 


DS  DETERMINES  THE  NO.  OF  POINTS  GENERATED  FOR  EACH 
TIME  HISTORY. 

DS=(UMAX-U(1))/TS 

NP  =  NO.  OF  FREQUENCIES  USED 


NP=1 


SPEC  98 
SPEC  99 
SPEC  100 
SPEC  101 
SPEC  102 
SPEC  103 
SPEC  104 
SPEC  105 
SPEC  106 
SPEC  107 
SPEC  108 
SPEC  109 
SPEC  110 
SPEC  111 
SPEC  112 
SPEC  113 
SPEC  114 
SPEC  115 
SPEC  116 
SPEC  117 
SPEC  118 
SPEC  119 
SPEC  120 
SPEC  121 
SPEC  122 
SPEC  123 
SPEC  124 
SPEC  125 
SPEC  126 
SPEC  127 
SPEC  128 
SPEC  129 
SPEC  130 
SPEC  131 
SPEC  132 
SPEC  133 
SPEC  134 
SPEC  135 
SPEC  136 
SPEC  137 
SPEC  138 
SPEC  139 
SPEC  140 
SPEC  141 
SPEC  142 
SPEC  143 
SPEC  144 
SPEC  145 
SPEC  146 
SPEC  147 
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C  FT  IS  FREQUENCY  TIMES  THE  PULSE  DURATION.  F  IS  FREQUENCY  SPEC 

C  AND  IS  EQUAL  TO  FT  SINCE  TIME  HAS  BEEN  NORMALIZED.  SPEC 

C  IF  IPLOT  EQUALS  1,  PLOT  TIME  HISTORY  OF  FREQUENCY  SPEC 

C  SPEC 

90  IF(IRES.EQ.1.0R.ISPEC.EQ.1.0R.I3D.EQ.l)  GO  TO  92  SPEC 

READ(5,1800)  FT, DAMP, IPLOT  SPEC 

WRITE(6,1800)  FT, DAMP, IPLOT  SPEC 

GO  TO  94  SPEC 

92  FT = FTR(NP)  SPEC 

IF(NP.EQ.l)  READ(5,1800)  DAMP  SPEC 

94  F0=FT  t,  .  -  ..  ..  ...  >  SPEC 

FT =FT*SQRT(1.-DAMP**2)  SPEC 

F= FT/UMAX  SPEC 

C  ,  SPEC 

C  XI  AND  X2  ARE  RUNNING  VALUES  OF  INTEGRATION  AND  ARE  SPEC 

C  SET  TO  0  WHEN  A  NEW  FREQUENCY  IS  READ.  SPEC 

C  SPEC 

XI =0.  SPEC 

X2=0.  SPEC 

1=1  SPEC 

K=1  SPEC' 

L=1  ,  SPEC 

M=1  SPEC 

c  }■»■■"...  : ;  ,  ..  SPEC 

C  H  IS  STEP  SIZE  FOR  SIMPSON’S  INTEGRATION  SPEC 

C  SPEC 

H=DS/8.  SPEC 

C  SPEC 

C  P  =  PERIOD  OF  FUNCTION.  SPEC 

C  SPEC 

P=l./F  SPEC 

C  SPEC 

C  H  MUST  BE  LESS  THAN  1/8  THE  PERIOD  FOR  REASONABLE  SPEC 

C  ACCURACY  USING  SIMPSON’S  INTEGRATION.  SPEC 

C  SPEC 

100  IF(H.LT.(1./8.)*P)  GO  TO  110  SPEC 

H  =  H/2.  SPEC 

GO  TO  100  SPEC 

110  H3=H/3.  SPEC 

C  SPEC 

C  SWS  =  SIN(WS)  AND  CWS  =  COS(WS).  SPEC 

C  SPEC 

W=PI2*F  SPEC 

WD=PI2*F0  SPEC 

120  WS=W*S  SPEC 

CWS  =  COS(WS)  SPEC 

SWS=SIN(WS)  SPEC 

C  SPEC 

C  U  AND  A  ARE  UPDATED  FOR  PROPER  LINEAR  INTERPOLATION.  SPEC 

C  SPEC 
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148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 

181 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 


130 

IF(T.GE.U(I)  AND.T.LE.U(I  + 1))  GO  TO  150 

SPEC  198 

IF(T.LT.UMAX)  GO  TO  140 

SPEC  199 

C 

SPEC  200 

c 

T  IS  SET  TO  TMAX  TO  CORRECT  POSSIBLE  ACCUMULATIVE  ERRORS 

SPEC  201 

c 

SPEC  202 

T=UMAX 

SPEC  203 

GO  TO  150 

SPEC  204 

140 

1=1  +  1 

SPEC  205 

GOTO  130 

SPEC  206 

c 

SPEC  207 

c 

Y  IS  THE  LINEAR  INTERPOLATION  F(T).  HENCE  THE  NEW 

SPEC  208 

c 

SET  OF  EQUALLY  SPACED  ORDERED  PAIRS  ARE  (T,Y) 

SPEC  209 

c 

SPEC  210 

150 

Y  =  (T-U(I))*(A(I  +  l)-A(I))/ (U(I  +  l)-U(I))  +  A(I) 

SPEC  211 

c 

SPEC  212 

c 

WT  =  2*PI*F*T  SWT  =  SIN(WT).  CWT  =  COS(WT). 

SPEC  213 

c 

SPEC  214 

WT  =W*T 

SPEC  215 

SWT=SIN(WT) 

SPEC  216 

CWT = COS(WT) 

SPEC  217 

c 

SPEC  218 

c 

THIS  SECTION  APPLIES  SIMPSON’S  INTEGRATION  FORMULA  TO 

SPEC  219 

c 

ORDERED  PAIRS  (T,V)  AND  (T,VP)  AND  STORES  THE  RESULTS 

SPEC  220 

c 

IN  X  AND  XP  RESPECTIVELY.  V  AND  VP  ARE  PARTS  OF 

SPEC  221 

c 

DUHAMEL’S  INTEGRAL  THAT  ARE  TO  BE  INTEGRATED. 

SPEC  222 

c 

SPEC  223 

V(K)=Y*SWT 

SPEC  224 

VP(K)=Y*CWT 

SPEC  225 

T =T+H 

SPEC  226 

IF(T.LT.0..OR.T.GT.UMAX + H)  GO  TO  300 

SPEC  227 

c 

SPEC  228 

c 

T  HAS  A  TOLERANCE  H/2  TO  ASSURE  THAT  LAST  POINT  IS  USED. 

SPEC  229 

c 

SPEC  230 

IF(T.GT.S  +  H/2.)  GO  TO  160 

SPEC  231 

K=K+1 

SPEC  232 

GO  TO  110 

SPEC  233 

c 

SPEC  234 

160 

OD=0. 

SPEC  235 

EV=0. 

SPEC  236 

DO  170  J=2,K,2 

SPEC  237 

OD=OD  + V(J) 

SPEC  238 

170 

EV=EV+ V(J  +  1) 

SPEC  239 

EV=EV-V(K) 

SPEC  240 

X = H3*(V(1)  +  V(K)+4.*OD + 2.*EV) 

SPEC  241 

OD=0. 

SPEC  242 

EV=0. 

SPEC  243 

DO  180  J=2,K,2 

SPEC  244 

OD=OD  + VP(J) 

SPEC  245 

180 

EV = E  V + VP(J  +  1) 

SPEC  246 

EV=EV-VP(K) 

SPEC  247 

38 


XP=H3*(VP(l)+VP*(K)+4.*OD+2.*EV) 

SPEC  248 

c 

SPEC  249 

c 

XI  AND  X2  SUM  TERMS  SO  THAT  WE  ARE  INTEGRATING  OVER 

SPEC  250 

c 

THE  RANGE  0.  TO  S. 

SPEC  251 

c 

SPEC  252 

X1=X1+X 

SPEC  253 

X2=X2+XP 

SPEC  254 

c 

SPEC  255 

c 

NX  IS  THE  INTEGRAL  OF  DUHAMEL’S  EQUATION  FOR  THE 

SPEC  256 

c 

RANGE  0.  TO  S. 

SPEC  257 

c 

SPEC  258 

NX  =  EXP(-D  AMP  *S  *  WD)  *W*(SWS*X2-CWS*X1) 

SPEC  259 

IF  (ABS(NX).LT..lE-06)  NX=0 

SPEC  260 

c 

i-i'r  ■■■:'.  i'  -hY  "  •  >t  •< 

SPEC  261 

c 

(T,NX)  IS  A  POINT  ON  THE  TIME  HISTORY  CURVE. 

SPEC  262 

c 

SPEC  263 

c 

WHEN  S  .LT.  UMAX  THE  PROGRAM  INCREASES  S  BY  DS  AND 

SPEC  264 

c 

DECREASES  T  BY  H.  IT  THEN  RETURNS  TO  LOCATION  120 

SPEC  265 

c 

AND  BEGINS  TO  INTEGRATE  THE  NEXT  AREA. 

SPEC  266 

c 

'  *  '  ’it.  » 

SPEC  267 

IF  (S.EQ.0.)  GO  TO  190 

SPEC  268 

1  =  1-1 

SPEC  269 

190 

T =T-H 

SPEC  270 

IF(IRES.EQ.l)  GO  TO  200 

SPEC  271 

c 

SPEC  272 

c 

TA  STORES  T 

SPEC  273 

c 

;  .u«  •"  '  V  '  • 

SPEC  274 

TA(L)  =  T 

SPEC  275 

c 

SPEC  276 

c 

XA  STORES  NX 

SPEC  277 

c 

SPEC  278 

XA(L)  =  NX 

SPEC  279 

XPLT(M)  =  T 

SPEC  280 

YPLT(M)  =  NX  ’ 

SPEC  281 

L=L+1  1 

SPEC  282 

M  =  M  +  1 

SPEC  283 

200 

K=1 

SPEC  284 

IF  ((S-UMAX).GE.-.1E-10)GO  TO  210 

SPEC  285 

S=S+DS 

SPEC  286 

GO  TO  120 

SPEC  287 

c 

SPEC  288 

c 

NXP  IS  THE  DERIVATIVE  OF  THE  TIME  HISTORY  AT  UMAX. 

SPEC  289 

c 

SPEC  290 

210 

ST  =  W*(T-S) 

SPEC  291 

SWST  =  SIN(ST) 

SPEC  292 

CWST  =  COS(ST) 

SPEC  293 

TT=T-S 

SPEC  294 

NXP  =  EXP(-DAMP*WD*TT)*W*W*(X1*SWST  +  X2*CWST) 

SPEC  295 

*  +EXP(-DAMP*WD*TT)*W*(DAMP*WD)*(-X2*SWST  +  X1*CWST) 

SPEC  296 

c 

SPEC  297 
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C  XMAX  IS  THE  MAXIMUM  VALUE  OF  THE  RESIDUAL  SPECTRUM 

C 

XMAX = SQRT(NX*NX + NXP*NXP/W/W) 

C 

L=1 

220  ST=WD*(T-S) 

SWST=SIN(ST) 

CWST = COS(ST) 

TT=T-S 

C 

C  (TB,XB)  ARE  ORDERED  PAIRS  FOR  THE  RESIDUAL  TIME  HISTORY 

C 

XB(L) = EXP(-DAMP*WD*TT)*(NXP*SWST/WD + NX*CWST) 

TB(L)=T 
XPLT(M)=T 
YPLT(M) =XB(L) 

M=M  +  1 
L=L+1 

IF(T.GE.3.*UMAX)  GO  TO  230 
T=T+DS 
GO  TO  220 
230  M=M-1 

IF(IPLOT.NE.l)  GO  TO  240 
READ(5,1100)  TIT4 
C 

C  PLOT  TIME  HISTORY  OF  INDIVIDUAL  FREQUENCY 

C 

LST=1 

LSP=M 

CALL  DATA4(NCURV,XPLT,YPLT,DUM,M,ISC,ITYPE,LST,LSP, 

*  TIT4,ICOLOR) 

C 

C  THIS  SECTION  FINDS  MAGNITUDE  OF  PRIMARY  TIME  HISTORY 

C  (MAXP)  AND  MAGNITUDE  OF  RESIDUAL  TIME  HISTORY  (MAXR) 

C 

240  MAXP=XA(1) 

MAXR=XB(1) 

DO  250  I  =  1,L 

IF(XA(I).LE.MAXP)  GO  TO  250 
MAXP=XA(I) 

250  CONTINUE 

DO  260  I  =  1,L 

IF(XB(I).LE.MAXR)  GO  TO  260 
MAXR=XB(I) 

260  CONTINUE 

C 

C  THIS  SECTION  STORES  DATA  FOR  3-D  PLOT 

C 

MSA  VE(NP) = M 

IF(NP.EQ.1AND.I3D.EQ.1)  WRITE(2)DS 
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SPEC  298 
SPEC  299 
SPEC  300 
SPEC  301 
SPEC  302 
SPEC  303 
SPEC  304 
SPEC  305 
SPEC  306 
SPEC  307 
SPEC  308 
SPEC  309 
SPEC  310 
SPEC  311 
SPEC  312 
SPEC  313 
SPEC  314 
SPEC  315 
SPEC  316 
SPEC  317 
SPEC  318 
SPEC  319 
SPEC  320 
SPEC  321 
SPEC  322 
SPEC  323 
SPEC  324 
SPEC  325 
SPEC  326 
SPEC  327 
SPEC  328 
SPEC  329 
SPEC  330 
SPEC  331 
SPEC  332 
SPEC  333 
SPEC  334 
SPEC  335 
SPEC  336 
SPEC  337 
SPEC  338 
SPEC  339 
SPEC  340 
SPEC  341 
SPEC  342 
SPEC  343 
SPEC  344 
SPEC  345 
SPEC  346 
SPEC  347 


IF  (I3D.EQ.1)  WRITE(2)  M,(XPLT(I),I  =  1,M),(YPLT(I),I  =  1,M) 

SPEC  348 

c 

SPEC  349 

IF(IFREQ.EQ.l)  F=(F*UMAX)/TCAL 

SPEC  350 

c 

SPEC  351 

c 

THIS  SECTION  CONVERTS  FREQUENCY  TIMES  PULSE 

SPEC  352 

c 

DURATION (FT)  TO  FREQUENCY,  IF  DESIRED,  AND  CONVERTS 

SPEC  353 

c 

FT  TO  LOG  10  SCALE,  IF  DESIRED,  AND  STORES  IT 

SPEC  354 

c 

ALONG  WITH  MAXR  AND  XMAX  FOR  PLOTTING. 

SPEC  355 

c 

'  r  •  ‘ 

SPEC  356 

IF(IFREQ.EQ.l)  FT=FT/TCAL 

SPEC  357 

IF(ILOG.NE.l)  PFT(NP)  =  FT 

SPEC  358 

IF(ILOG.EQ.l)  PFT (NP)  = ALOG  10(FT) 

SPEC  359 

RMAX(NP)  =  MAXR 

SPEC  360 

FR(NP)  =  F 

SPEC  361 

FWN(NP)  =  FO/TCAL 

SPEC  362 

PXP(NP)  =  MAXP 

SPEC  363 

PXR(NP) = MAXR 

SPEC  364 

NP=NP+1 

SPEC  365 

c 

SPEC  366 

c 

AFTER  THE  CODE  COMPUTES  THE  SPECTRUM  FOR  FMAX,  IT 

SPEC  367 

c 

GOES  TO  THE  PLOT  SUBROUTINE.  OTHERWISE,  IT  RESETS  T, 

SPEC  368 

c 

I,  AND  S  AND  RETURNS  TO  LOCATION  90. 

SPEC  369 

c 

SPEC  370 

IF(FO.GEJFMAX)  GO  TO  270 

SPEC  371  - 

T =0 

SPEC  372 

1  =  1 

SPEC  373 

s=o 

SPEC  374 

GO  TO  90 

SPEC  375 

c 

SPEC  376 

c 

NP  IS  NO.  OF  FREQUENCIES  USED.  IT  IS  RESTRICTED  TO 

SPEC  377 

c 

200  DUE  TO  DIMENSION  STATEMENT 

SPEC  378 

c 

-  * . 

SPEC  379 

270 

NP=NP-1 

SPEC  380 

c 

SPEC  381 

c 

PLOT  PRIMARY  AND  RESIDUAL  SPECTRA 

SPEC  382 

c 

SPEC  383 

IF(ISPEC.NE.l)  GO  TO  280 

SPEC  384 

READ(5,1100)  TIT4 

SPEC  385 

ITYPE=1 

SPEC  386 

LST=1 

SPEC  387 

LSP  =  NP 

SPEC  388 

NCURV=2 

SPEC  389 

ISC=1 

SPEC  390 

IF(ISPEC.EQ.1)CALL  DATA4(NCURV,PFT,PXP,PXR,NP,ISC, 

SPEC  391 

*ITYPE,LST,LSP,TIT4,ICOLOR) 

SPEC  392 

IF(IRES.EQ.l)  NCURV=1 

SPEC  393 

IF(IRES.EQ.l)  CALL  DATA4(NCURV,PFT,PXR,DUM,NP,ISC, 

SPEC  394 

*ITYPE,LST,LSP,TIT4,ICOLOR) 

SPEC  395 

ITYPE=0 

SPEC  396 

ISC=0 

SPEC  397 

41 


c 

SPEC  398 

c 

PRINT  HEADINGS  AND  DATA 

SPEC  399 

c 

SPEC  400 

280 

WRITE(6,1200)  PGTIT 

SPEC  401 

WRITE(6,1900) 

SPEC  402 

WRITE(6,2000) 

SPEC  403 

IF(IFREQ.EQ.l)  WRITE(6,2100) 

SPEC  404 

IF(IFREQ.NE.l)  WRITE(6,2200) 

SPEC  405 

WRITE(6,2500)  (FR(I),PXP(I),RMAX(I),FWN(I),I= 1,NP) 

SPEC  406 

WRITE(6,2300)  NP 

SPEC  407 

REWIND  2 

SPEC  408 

DO  290  1  =  1,1000 

SPEC  409 

V(I)=0.0 

SPEC  410 

VP(I)=0. 

SPEC  411 

290 

CONTINUE 

SPEC  412 

GO  TO  10 

SPEC  413 

300 

WRITE(6,2400)  T 

SPEC  414 

STOP 

SPEC  415 

1000 

FORMAT  (2I5,3F10.3,F5.0,5I5,3I2) 

SPEC  416 

1100  . 

FORMAT(10A8) 

SPEC  417 

1200 

FORMAT(1H1,10A8,/) 

SPEC  418 

1300 

FORMAT(1HO,4(6X, TIME’, 7X, ’ACCELERATION’, IX)) 

SPEC  419 

1400 

FORMAT(  1H  ,4(6X,4HMSEC,1 1X,3HG’S,6X) /) 

SPEC  420 

1500 

FORMAT  (8(1X,F14.5)) 

SPEC  421 

1600 

FORMAT(8F10.3) 

SPEC  422 

1700 

FORMAT(lHl, ’NORMALIZED  INPUT  DATA’/) 

SPEC  423 

1800 

FORMAT(2F10.3,I5) 

SPEC  424 

1900 

FORMAT(’  PRIMARY  AND  RESIDUAL  SPECTRA  OF 

SPEC  425 

♦ACCELERATION  PULSE’,//) 

SPEC  426 

2000 

FORMAT(’  FREQUENCY’, 5X,’  RELATIVE 

SPEC  427 

♦  ACCELERATION  ACTUAL’) 

SPEC  428 

2100 

FORMAT(’  KHZ’, 13X, ’PRIMARY  RESIDUAL 

SPEC  429 

♦FREQUENCY’) 

SPEC  430 

2200 

FORMAT(’  X  PULSE  DURATION’, 5X,’  PRIMARY 

SPEC  431 

♦  RESIDUAL  FREQUENCY’) 

SPEC  432 

2300 

FORMAT(’0NFILES  =  ’,15) 

SPEC  433 

2400 

FORMAT(’0T<0  OR  >  UMAX  +  H,  T=’,F10.3) 

SPEC  434 

2500 

FORMAT  (F12.3,5X,3F12.5) 

SPEC  435 

C 

END 

SPEC  436 

c 

c 


SUBROUTINE  DATA4(NCURV,Z1,Z2,Z3,NPLT, ISC, ITYPE.LST, 

DAT4  1 

♦LSP, TITLE, ICOLOR) 

DAT4  2 

C 

DAT4  3 

C 

THIS  SUBROUTINE  SETS  UP  DATA  AND  LABELS  FOR  PLOT 

DAT4  4 

C 

DAT4  5 

DIMENSION  Z1(NPLT),Z2(NPLT),X(1000),Y1(1000),DUM(1), 

DAT4  6 

*Y2(1000),Z3(NPLT),LXNAME(3),LYNAME(3),TITLE(4) 

DAT4  7 

DUM(1)  =  1.E100 

DAT4  8 

42 


o  n 


c 

C  STORE  DATA  TO  BE  PLOTTED  IN  APPROPRIATE  ARRAYS  AND 

C  COUNT  DATA  POINTS 

C 

J=0 

DO  100  L=LST,LSP 
J=J  +  1 
X(J)  =  Z1(L) 

Y1(J)  =  Z2(L) 

Y2(J)=Z3(L) 

100  CONTINUE 

IF(NCURV.LT.2)  Y2(1)=DUM(1) 

C 

C  PRINT  PLOT  LIMITS 

C 

WRITE(6,2)  LST,LSP,X(1),X(J) 

C 

C  READ  AND  PRINT  PLOT  LABELS 

C 

READ(5,3)  LXNAME.LYNAME 
WRITE(6,3)  LXNAME.LYNAME 
C 

C  CALL  PLOT  SUBROUTINE 

C 

CALL  DISS4(X,Y1,Y2,J, TITLE, ISC, ITYPE,LXNAME,LYNAME, 
*ICOLOR) 

RETURN 

1  FORMAT  (2F10.0) 

2  FORMAT(’  PLOT  LIMITS  -  \2I5,2F10.4) 

3  FORMAT(10A8) 

END 


SUBROUTINE  DISS4(X,Yl,Y2,NPT, TITLE, ISC, ITYPE.LXNAME, 
*LYNAME,ICOLOR) 

C 

C  THIS  SUBROUTINE  CREATES  A  DISSPLA  9  PLOT  FILE 

C 

C 

C  ISC=0  SELF-SCALE 

C  =1  READ  IN  XORIG,  XMAX,  XSTP,  YORIG,  YMAX,  YSTP 

C 

C  ITYPE=0  NORMAL  PLOT 

C  =  1  X-LOG  AXIS 

C 

C  ICOLOR=0  BLACK/WHITE 

C  =1  COLOR 

C 

DIMENSION  X(NPT),Y1(NPT),Y2(NPT),RAT(10),LXNAME(3), 

*  LYNAME(3),TITLE(4) 


DAT4  9 
DAT4  10 
DAT4  11 
DAT4  12 
DAT4  13 
DAT4  14 
DAT4  15 
DAT4  16 
DAT4  17 
DAT4  18 
DAT4  19 
DAT4  20 
DAT4  21 
DAT4  22 
DAT4  23 
DAT4  24 
DAT4  25 
DAT4  26 
DAT4  27 
DAT4  28 
DAT4  29 
DAT4  30 
DAT4  31 
DAT4  32. 
DAT4  33 
DAT4  34 
DAT4  35  - 
DAT4  36 
DAT4  37 
DAT4  38 
DAT4  39 


DIS4 

1 

DIS4 

2 

DIS4 

3 

DIS4 

4 

DIS4 

5 

DIS4 

6 

DIS4 

7 

DIS4 

8 

DIS4 

9 

DIS4 

10 

DIS4 

11 

DIS4 

12 

DIS4 

13 

DIS4 

14 

DIS4 

15 

DIS4 

16 

DIS4 

17 

43 


c 

DIS4  18 

c 

DETERMINE  PLOT  SCALES 

DIS4  19 

c 

DIS4  20 

IF(ISC.NE.O)GO  TO  300 

DIS4  21 

XORIG=X(l) 

DIS4  22 

XSTP= ’SCALE’ 

DIS4  23 

XMAX = X(NPT) 

DIS4  24 

YORIG  =  1.E100 

DIS4  25 

YSTP= ’SCALE’ 

DIS4  26 

YMAX=-1.E100 

DIS4  27 

DO  200  I  =  1,NPT 

DIS4  28 

IF(Yl(I).LT.YMAX)GO  TO  100 

DIS4  29 

YMAX=Y1(I) 

DIS4  30 

100 

IF(Y  1(I).GT.Y ORIG)GO  TO  200 

DIS4  31 

YORIG=Yl(I) 

DIS4  32 

200 

CONTINUE 

DIS4  33 

GO  TO  400 

DIS4  34 

300 

READ(5,l)XORIG,XMAX,XSTP,YORIG,YMAX,YSTP 

DIS4  35 

400 

IF(ISC.EQ.l)WRITE(6,2)XORIG,XSTP,XMAX,YORIG,YSTP,YMAX 

DIS4  36 

IF(ISC.EQ.O)WRITE(6,2)XORIG,XMAX,YORIG,YMAX 

DIS4  37 

JJ=JJ  +  1 

DIS4  38 

c 

DIS4  39 

c 

SET  UP  GRID 

DIS4  40 

c 

DIS4  41 

CALL  BASALF(’STANDARD’) 

DIS4  42 

CALL  MIXALF(’SPECIAL’) 

DIS4  43 

CALL  PHYSOR(1.5,l.) 

DIS4  44 

CALL  PAGE(U.,8.5) 

DIS4  45 

CALL  SETCLR(’BLACK’) 

DIS4  46 

XSIZE=8.0 

DIS4  47 

YSIZE=6.0 

DIS4  48 

CALL  AREA2D(XSIZE,YSIZE) 

DIS4  49 

c 

DIS4  50 

c 

LABEL  PLOTS 

DIS4  51 

c 

DIS4  52 

CALL  HEIGHT(.25) 

DIS4  53 

CALL  HEADIN(TITLE, 32, 1.1,1) 

DIS4  54 

CALL  XNAME(LXNAME,24) 

DIS4  55 

CALL  YNAME(LYNAME,24) 

DIS4  56 

c 

DIS4  57 

c 

SET  PLOT  SCALES  FOR  LOG  PLOT 

DIS4  58 

c 

DIS4  59 

IF(ITYPE.NE.l)  GO  TO  500 

DIS4  60 

ICY = ALOG10(X(NPT))-ALOG10(X(1))  +  .5 

DIS4  61 

XCY  =  XSIZE/FLOAT(ICY) 

DIS4  62 

IF(ITYPE.EQ.1AND.ISC.EQ.O)YSTP=(YMAX-YORIG)/YSIZE 

DIS4  63 

WRITE(6,3)  XCY,YSTP,ICY 

DIS4  64 

c 

D1S4  65 

c 

DRAW  GRID 

DIS4  66 

c 

DIS4  67 

44 


n  n  n 


500 


IF(ITYPE.EQ.O)CALL  GRAF(XORIG,XSTP,XMAX,YORIG, 

DIS4 

68 

♦YSTP.YMAX) 

DIS4 

69 

IF(ITYPE.EQ.l)CALLXLOG(XORIG^CCY,YORIG,YSTP) 

DIS4 

70 

IMARK=0 

DIS4 

71 

IF(ICOLOR.EQ.l)CALL  SETCLR(’BLUE’) 

DIS4 

72 

DIS4 

73 

PLOT  DATA 

DIS4 

74 

DIS4 

75 

CALL  CURVE(X,Y1,NPT,IMARK) 

DIS4 

76 

IF(ICOLOR.EQ.l)CALL  SETCLR(’RED’) 

DIS4 

77 

TL=.5 

DIS4 

78 

NMRK=10 

DIS4 

79 

RAT(1)=2. 

DIS4 

80 

RAT(2)  =  1. 

DIS4 

81 

IF(ICOLOR.EQ.O)CALL  MRSCOD(TL,NMRK,RAT) 

DIS4 

82 

IF(Y2(1).NE.1.E100AND.ITYPE.EQ.0) 

DIS4 

83 

*CALL  CURVE(X,Y2,NPT,IMARK) 

DIS4 

84 

IF(Y2(l).NE.l.E100AND.ITYPE.EQ.l) 

DIS4 

85 

*CALL  CURVE(X,Y2,NPT,IMARK) 

DIS4 

86 

CALL  RESET(’DASH’) 

DIS4 

87 

CALL  ENDPL(JJ) 

DIS4 

88 

RETURN 

DIS4 

89 

FORMAT (6E10.3) 

DIS4 

90 

FORMAT(’  PLOT  SCALES  -  ’.6E13.5) 

DIS4 

91- 

FORMAT(’  XCYCLE,  YSTP,  ICY  = ’,2E13.5,I10) 

DIS4 

92 

END 

DIS4 

93 
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APPENDIX  B 


LISTING  OF  FORTRAN  PROGRAM  TO  CREATE  3-D  SHOCK  RESPONSE  PLOT 
ALONG  WITH  ASSOCIATED  JOB  CONTROL  LANGUAGE  FOR  CRAY  XMP/48 


I.  Job  Control  Language  for  CRAY  XMP/48 


JOB,JN = jobname  ,T  =  150,MFL= 350000. 

ACCOUNT  AC = accountnumber  ,US = username  ,UPW =userpassword. 

ACCESS, DN = DISLIB,PDN = DISLIB.ID = DISSPLA,OWN= SYSTEM. 

ACCESS, DN = DVSD,PDN = DVSD,ID = DISSPLA,0  WN = SYSTEM. 

ACCESS, DN = INTLIB,PDN = INTLIB.ID = DISSPLA.O  WN = SYSTEM. 

FETCH, DN=3dfile; TEXT =' CHARGE ficctno ,  pn  .GET,  3dfile  =3dfile  ,CTASK.’,DF=BB. 

where  3 dfile  =  name  of  input  file  which  contains  data  for 
calculations  for  3-D  plot 
pn  —  first  two  letters  of  account  number 

ASSIGN, DN=3dfileA=FT01. 

FETCH, DN=/n/i/e  ,TEXT = ’CHARGE pcctno ,  pn  .GET,  infile  =infile  .CTASK.’. 

where  infile  =  name  of  input  file  which  contains  program  options 
and  plot  titles 

pn  =  first  two  letters  of  the  account  number 

ASSIGN, DN=infile  ,A=FT05. 

CFT. 

SEGLDR.CMD  =  ’LIB  =  DISLIB,INTLIB’,GO. 

DISPOSE, DN = META, DC = ST,DF=BB,/S 

TEXT = ’CHARGE = acctno ,  pn  ,CTASK.REPLACE,META  =pl  3 file .’. 

where  pi  3 file  =  name  of  output  file  in  which  the  DISSPLA  plot  file 
for  3-D  plot  is  stored  on  the  front  end. 


II.  Listing  of  FORTRAN  File 


PROGRAM  DUHM3D(INPUT, OUTPUT, TAPE1, TAP E5,  DU3D  1 

*  TAP  E6 = OUTPUT, MET  A, TAPE10 = OUTPUT, TAPE2)  DU3D  2 

C  DU3D  3 

C  THIS  PROGRAM  READS  FREQUENCY  RESPONSE  DATA  AND  DU3D  4 

C  CREATES  A  3-D  SURFACE  PLOT  DU3D  5 

C  DU3D  6 

C  THIS  PROGRAM  USES  A  DATA  FILE  CREATED  BY  A  DU3D  7 

C  DUHAMEL  INTEGRAL  PROGRAM  AS  INPUT  DU3D  8 
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C  DU3D 

DIMENSION  TX(65),TY (65),TX1(65),TY  1(65)  DU3D 

DIMENSION  TX2(65),TY2(65),TX3(65),TY3(65),TX4(65),TY4(65)  DU3D 

DIMENSION  OFF(50),OFFX(50),PARRAY(1500,4),RES(1500,4),  DU3D 

*  T1(4,4),XPT(1500),YPT(1500),TIT4(4),XPLT(1500),YPLT(1500),  DU3D 

*  T2(4,4)  DU3D 

C  DU3D 

DATA  OFF/L,L33,L67, 2,233.3.5, 4,5,6,7,8,9,10,13.3,  DU3D 

♦16.7, 20., 25., 30., 35., 40,50., 60., 70..80..90., 100., 13333,166.67,  DU3D 

*200,250,300,350,400,500,600,700,800,900,1000,1333.3,  DU3D 

*1666.67,2000,2500,3000,4000,5000,6000,7000,8000./  DU3D 

DU3D 

INITIALIZE  DISSPLA  DU3D 

DU3D 

CALL  COMPRS  DU3D 

CALL  SETDEV(10,6)  DU3D 

C  DU3D 

ISC=1.  DU3D 

ICOLOR  =  1  DU3D 

IEP=1  DU3D 

C  DU3D 

C  NFILES(I5)  =  NUMBER  OF  FREQUENCY  RESPONSE  CURVES  DU3D 

C  MAXT  (15)  =  MAXIMUM  TIME  OF  RESPONSE  PLOT  DU3D 

C  MAXT  =  1,  PRIMARY  RESPONSE  CURVE  DU3D 

C  MAXT  =  2,  PRIMARY  RESPONSE  PLUS  PARTIAL  RESIDUAL  DU3D 

C  MAXT  =  3,  PRIMARY  RESPONSE  PLUS  FULL  RESIDUAL  DU3D 

C  DU3D 

READ(5,1000)  NFILES.MAXT  DU3D 

NF=NFILES  DU3D 

C  DU3D 

C  DS  -  TIME  STEP  READ  FROM  DATA  FILE  DU3D 

C  DU3D 

READ(l)  DS  DU3D 

WRITE(6,5000)  DS  DU3D 

C  DU3D 

C  CALCULATE  NUMBER  AND  INTERVAL  OF  TIME  CURVES  DU3D 

C  TO  BE  PLOTTED  DU3D 

C  DU3D 

K = (MAXT/DS)  + 1  DU3D 

A  =  MAXT  DU3D 

B=AMOD(A,DS)  DU3D 

IF(B.GE..0)K=K+1  DU3D 

KD  =  K/4  DU3D 

Jl=4  DU3D 

J  =  MOD(KrJl)  DU3D 

IF(J.GT.0)KD  =  KD  + 1  DU3D 

KS=K/100  DU3D 

Jl=100  DU3D 

J  =  MOD(K31)  DU3D 

IF(J.GT.0)KS = KS  + 1  DU3D 


9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 
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KD2=2*KD 

DU3D  59 

KD3=3*KD 

DU3D  60 

KP=0 

DU3D  61 

WRITE(6,6000)  KD,KS,K 

DU3D  62 

IPSAVE  =  10000 

DU3D  63 

c 

DU3D  64 

c 

READ  DATA  AND  STORE  IN  SEPARATE  ARRAYS  TO  CREATE 

DU3D  65 

c 

DATA  CURVES  AT  SPECIFIC  TIME  INTERVALS 

DU3D  66 

c 

DU3D  67 

DO  80  J  =  1,KD,KS 

DU3D  68 

IS=0 

DU3D  69 

DO  60  1  =  1, NF 

DU3D  70 

IS=IS  +  1 

DU3D  71 

IF(J.GT.lAND.I.EQ.l)  READ(l)  DS 

DU3D  72 

READ(l)  M,(XPLT(N),N = 1,M),(YPLT(N),N = 1,M) 

DU3D  73 

IF(I.EQ.l  AND  J.EQ.1)WRITE(6, 4000)  (XPLT(N),N= 1,M),(YPLT(N),N= 1,M) 

DU3D  74 

IF(J.EQ.1AND.I.EQ.1)  GO  TO  30 

DU3D  75 

10 

IF((J + KD2)  .GT.IPS A VE)  KP  =  1 

DU3D  76 

TX(I)  =XPLT(J) 

DU3D  77 

TY  (I)  =  YPLT(J) 

DU3D  78 

TX1(I)=XPLT(J +KD) 

DU3D  79 

TY1(I)  =  YPLT(J +KD) 

DU3D  80 

TX2(I) = XPLT(J  +  KD2 + KP) 

DU3D  81 

TY2(I) = YPLT(J + KD2 + KP) 

DU3D  82 

TX3(I) = XPLT(J + KD3 + KP) 

DU3D  83 

TY3(I) = YPLT(J + KD3 + KP) 

DU3D  84 

IF((J  +  KS).GT.KD AND  J.NE.KD)  GO  TO  20 

DU3D  85 

GO  TO  60 

DU3D  86 

20 

IF(I.EQ.l)  IEND  =  1 

DU3D  87 

TX4(I)  =XPLT(K) 

DU3D  88 

TY4(I)  =  YPLT(K) 

DU3D  89 

GO  TO  60 

DU3D  90 

30 

Kl=l 

DU3D  91 

DO  40  IP  =  1,M 

DU3D  92 

IF((XPLT(IP  +  1)-XPLT(IP)).LT..001)  IPSAVE= IP 

DU3D  93 

IF(ABS(MAXT-XPLT(IP)).LT..0001)  GO  TO  50 

DU3D  94 

K1=K1+1 

DU3D  95 

40 

CONTINUE 

DU3D  96 

50 

K=K1 

DU3D  97 

GO  TO  10 

DU3D  98 

60 

CONTINUE 

DU3D  99 

REWIND  1 

DU3D  100 

C 

DU3D  101 

c 

STORE  NEWLY  CREATED  DATA  CURVES 

DU3D  102 

c 

DU3D  103 

WRITE(2)  IS,(TX(N),N = 1,IS),(TY(N),N = 1,IS) 

DU3D  104 

NFILES = NFILES  + 1 

DU3D  105 

WRITE(2)  (TX1(N),N = 1,IS),(TY1(N),N = 1,IS) 

DU3D  106 

NFILES = NFILES  + 1 

DU3D  107 

WRITE (2)  (TX2(N),N = 1,IS),(TY2(N),N = 1,IS) 

51 

DU3D  108 

u  u  u 


NF1LES  =  N  FILES + 1  DU3D 

WRITE(2)  (TX3(N),N = 1,IS),(TY3(N),N = 1,IS)  DU3D 

NFILES = NFILES  + 1  DU3D 

WRITE(6,4000)  (TX(N),N  =  1,8),(TY(N),N = 1,8)  DU3D 

IF  (IEND.NE.1)  GO  TO  70  DU3D 

WRITE(2)  (TX4(N),N = 1,IS),(TY4(N),N = 1,IS)  DU3D 

NFILES = NFILES + 1  DU3D 

70  WRITE(6,2000)  J  DU3D 

IF  (IEND.EQ.l)  WRITE(6,3000)  DU3D 

80  CONTINUE  DU3D 

REWIND  2  DU3D 

DU3D 

GENERATE  OFFSET  FACTORS  FOR  PLOTTING  DU3D 

DU3D 

RMS  =  (ALOG10(OFF(NF))-ALOG10(OFF(1)))/(NF*.01)  DU3D 

DO  90  1  =  1,NF  DU3D 

OFFX(I)  =  (ALOG10(OFF(NF))-ALOG10(OFF(I))-NF*.01*RMS)/(-RMS)  DU3D 

90  CONTINUE  DU3D 

C  DU3D 

C  SET  UP  ARRAYS  FOR  CURVE  ROTATION  DU3D 

C  DU3D 

AX=1HY  DU3D 

ANG=53  DU3D 

CALL  ROT1  (ANG,AX,T1)  DU3D 

AX=1HX  DU3D 

ANG=67  DU3D 

CALL  ROT1  (ANGAX.T2)  DU3D 

C  DU3D 

C  READ  DATA  FROM  FILES  DU3D 

C  DU3D 

IMAXC  =  NFILES  DU3D 

DO  180  II  =  l.NFILES  DU3D 

Jl=4  DU3D 

J  =  MOD((II-NF),Jl)  DU3D 

MM=1  DU3D 

IF(II.GT.NF)  GO  TO  100  DU3D 

IF(II.EQ.l)  READ(l)  DS  DU3D 

READ(l)  M,(XPLT(I),I  =  1,M),(YPLT(I),I  =  1,M)  DU3D 

GO  TO  120  DU3D 

100  IF(J.EQ.lAND.II.NE.NFILES)GO  TO  110  DU3D 

READ(2)  (XPLT(I),I  =  1,M),(YPLT(I),I = 1,M)  DU3D 

GO  TO  120  DU3D 

110  READ(2)  M,(XPLT(I),I  =  1,M),(YPLT(I),I  =  1,M)  DU3D 

C  DU3D 

C  ROTATE  DATA  DU3D 

C  DU3D 

120  DO  130  JJ  =  1,M  DU3D 

PARRAY  (MM,2) = YPLT(JJ)  DU3D 

PARRAY  (MM,1) = XPLT(JJ)  DU3D 

PARRAY  (MM,3) = 0.  DU3D 


109 

110 
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129 
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147 
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150 
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158 


52 


PARRAY  (MM, 4)  =  1. 

DU3D  159 

MM  =  MM  + 1 

DU3D  160 

IF(II.GT.NF)  GO  TO  130 

DU3D  161 

IF(XPLT(JJ).GE.MAXT)  GO  TO  140 

DU3D  162 

130 

CONTINUE 

DU3D  163 

140 

MM  =  MM-1 

DU3D  164 

IA=MM 

DU3D  165 

IM  =  1500 

DU3D  166 

JM=4 

DU3D  167 

KM=4 

DU3D  168 

NM= 1500 

DU3D  169 

CALLMATMPY(PARRAY,T1,RES,IA,JM)KM,IM,JM,NM) 

DU3D  170 

CALL  MATMPY (RES, T2, PARRAY, IArFM,KM,IMrIM,NM) 

DU3D  171 

C 

DU3D  172 

c 

OFFSET  DATA  FOR  PLOTTING 

DU3D  173 

c 

DU3D  174 

DO  170  JJ  =  1,MM 

DU3D  175 

IF(II.LE.NF)  GO  TO  150 

DU3D  176 

OFFSETX = OFFX(JJ) 

DU3D  177 

OFFSET =OFF(JJ) 

DU3D  178 

IF(IEND.EQ.l  AND.II.EQ.NFILES  AND  JJ.EQ.MM)OFFSET = OFF(50) 

DU3D  179 

GO  TO  160 

DU3D  180 

150 

OFFSETX =OFFX(II) 

DU3D  181 

OFFSET =OFF(II) 

DU3D  182 

160 

XPT(JJ)  =  PARRAY(JJ.l)  +  OFFSETX 

DU3D  183 

YPT(JJ)  =  PARRAY  (JJ,2)-ALOG10(OFFSET) 

DU3D  184 

170 

CONTINUE 

DU3D  185 

C 

DU3D  186 

c 

READ  PLOT  TITLE  AND  CALL  PLOT  ROUTINE 

DU3D  187 

c 

DU3D  188 

IF(IEP.EQ.l)  READ(5,7000)  TIT4 

DU3D  189 

LST=1 

DU3D  190 

LSP=MM 

DU3D  191 

■  * 

NPLT = 1500 

DU3D  192 

CALL  DATA4(XPT)YPT)NPLT)ISC,LST,LSP,TIT4,ICOLOR,IEP,IMAXC) 

DU3D  193 

IEP=IEP+1 

DU3D  194 

180 

CONTINUE 

DU3D  195 

CALL  DONEPL 

DU3D  196 

1000 

FORMAT(2I5) 

DU3D  197 

2000 

FORMAT(’  CYCLE  NO.  ’,15,’  HAS  BEEN  STORED’) 

DU3D  198 

3000 

FORMAT(’  FINAL  FILE  HAS  BEEN  STORED’) 

DU3D  199 

4000 

FO  RMAT  (8F12.5) 

DU3D  200 

5000 

FORMAT(’0DS  =  \F30.20) 

DU3D  201 

6000 

FORMAT(’  KD  KS  K  ’,3110) 

DU3D  202 

7000 

FORMAT  ( 10A8) 

DU3D  203 

STOP 

DU3D  204 

END 

DU3D  205 

c 

SUBROUTINE  ROTl(ANG,AX,T) 

ROT1  1 
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c 

ROT1  2 

c 

THIS  SUBROUTINE  CREATES  A  MATRIX  FOR  ROTATION  OF  DATA 

ROT1  3 

c 

ROT1  4 

c 

ROT1  5 

DIMENSION  T(4,4) 

ROT1  6 

c 

ROT1  7 

c 

CONVERT  DEGREES  TO  RADIANS 

ROT1  8 

c 

ROT1  9 

ANG  =ANG*. 01745329 

ROT1  10 

c 

ROT1  11 

c 

ZEROING  OUT  THE  T  ARRAY 

ROT1  12 

c 

ROT1  13 

DO  1001=1,4 

ROT1  14 

DO  100  J= 1,4 

ROT1  15 

T(I,J)=0.0 

ROT1  16 

100 

CONTINUE 

ROT1  17 

IF(AX.EQ.1HY)  GO  TO  200 

ROT1  18 

IF(AX.EQ.1HZ)  GO  TO  300 

ROT1  19 

c 

ROT1  20 

c 

X-AXIS  ARRAY 

ROT1  21 

c 

ROT1  22 

T(2,3) = SIN(ANG) 

ROT1  23 

T  (2,2) = COS(ANG) 

ROT1  24 

T(3,2)  =  -T(2,3) 

ROT1  25 

T(3,3)=T(2,2) 

ROT1  26 

T(l,l)=l. 

ROT1  27 

GO  TO  400 

ROT1  28 

c 

ROT1  29 

c 

Y-AXIS  ARRAY 

ROT1  30 

c 

ROT1  31 

200 

T(3,l) = SIN(ANG) 

ROT1  32 

T(l,l) =COS(ANG) 

ROT1  33 

T(l,3)=-T(3,l) 

ROT1  34 

T(3,3)=T(1,1) 

ROT1  35 

T(2,2)  =  l. 

ROT1  36 

GO  TO  400 

ROT1  37 

c 

ROT1  38 

c 

Z-AXIS  ARRAY 

ROT1  39 

c 

ROT1  40 

300 

T(l,3) =SIN(ANG) 

ROT1  41 

T(l,l)=COS(ANG) 

ROT1  42 

T(2,l)  =  -T(l,3) 

ROT1  43 

T(2,2)=T(1,1) 

ROT1  44 

T(3,3)  =  l. 

ROT1  45 

400 

RETURN 

ROT1  46 

END 

ROT1  47 

c 

SUBROUTINE  DATA4(Z1,Z2,NPLT, ISC, LST,LSP, TITLE, ICOLOR, 

DAT4  1 

*IEP,IMAXC) 

DAT4  2 
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n  n  n  5  rro  o  n  o  o  o  o.o.p  n  n  o  n  n 


c 

c 

c 


20 

30 

40 

50 


C 


THIS  SUBROUTINE  SETS  UP  DATA  FOR  PLOTTING 
AND  CREATES  A  3-D  DISSPLA  9  PLOT  FILE 

DIMENSION  Z1(NPLT),Z2(NPLT) ,X(1000), Y 1(1000) ,TITLE(4) 

ISC  =  0  SELF-SCALE 

=  1  READ  IN  XORIG,  XMAX,  XSTP,  YORIG,  YMAX,  YSTP 

ICOLOR=0  BLACK/WHITE 
=  1  COLOR 


STORE  DATA  TO  BE  PLOTTED  IN  APPROPRIATE  ARRAYS  AND 
COUNT  DATA  POINTS 

'■ . i >?■■,  ..  £  ;»  .*«  ■ 

j=°  .  ? •  >.  ,  . 

DO  10  L=LST,LSP 

J=J  +  1 
X(J)=Z1(L) 

Y1(J)  =  Z2(L) 

CONTINUE 

PRINT  PLOT  LIMITS 

WRITE(6,3000)LST,LSP,X(1),X(J) 

NPT =J 

IF(IEP.GT.l)  GO  TO  60 

DETERMINE  PLOT  SCALES 

IF(ISC.NE.O) GO  TO  40 
XORIG =X(1) 

XSTP = ’SCALE’ 

XMAX=X(NPT) 

YORIG  =  1.E100 
YSTP  =  ’SCALE’ 

YMAX  =  -1.E100 
DO  30  1=  1,NPT 
IF(Yl(I).LT.YMAX)GO  TO  20 
YMAX=Y1(I) 

IF(Yl(I).GT.YORIG)GO  TO  30 
YORIG  =  Y1(I) 

CONTINUE 
GO  TO  50 

READ(5,1000)XORIG,XMAX, XSTP, YORIG, YMAX, YSTP 
IF(ISC.EQ.1)WRITE(6,2000)XORIG, XSTP, XMAX, YORIG, YSTP, YMAX 
IF(ISC.EQ.0)WRITE(6,2000)XORIG, XMAX, YORIG, YMAX 
JJ=JJ  +  1 


DAT4  3 
DAT4  4 
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c 

SET  UP  PLOTTING  AREA 

DAT4  53 

c 

DAT4  54 

CALL  PHYSOR(1.5,l.) 

DAT4  55 

CALL  NOBRDR 

DAT4  56 

CALL  PAGE(11.,8.5) 

DAT4  57 

CALL  SETCLR(’BLACK’) 

DAT4  58 

XSIZE=9.0 

DAT4  59 

YSIZE=6.0 

DAT4  60 

CALL  AREA2D(XSIZE,YSIZE) 

DAT4  61 

c 

DAT4  62 

c 

PLOT  TITLE 

DAT4  63 

c 

DAT4  64 

CALL  HEIGHT(.25) 

DAT4  65 

CALL  HEADIN(TITLE, 32, 1.1,1) 

DAT4  66 

CALL  GRAF(XORIG,XSTP,XMAX,YORIG,YSTP,YMAX) 

DAT4  67 

IMARK=0 

DAT4  68 

IF(ICOLOR.EQ.l)CALL  SETCLR(’BLUE’) 

DAT4  69 

c 

DAT4  70 

c 

PLOT  DATA 

DAT4  71 

c 

DAT4  72 

60 

CALL  CUR VE(X,Y  1,NPT,IMARK) 

DAT4  73 

IF(IMAXC.EQ.IEP)  CALL  ENDPL(JJ) 

DAT4  74 

RETURN 

DAT4  75 

1000 

FORMAT  (6E10.3) 

DAT4  76 

2000 

FORMAT(’  PLOT  SCALES  -  ’.6E13.5) 

DAT4  77 

3000 

FORMAT(’  PLOT  LIMITS  -  ’.2I5.2F10.4) 

DAT4  78 

r 

END 

DAT4  79 

V-/ 

SUBROUTINE  MATMPY  (A,B,C,U,K,L,M,N) 

MATMPY 

c 

MATMPY 

c 

A(U)  B(J,K)  ARE  ACTUAL  DIMENSIONS  WITHIN 

MATMPY 

c 

MAX  DIMENSIONS  OF  A(L,)  B(M,)  C(N,) 

MATMPY 

c 

REST  OF  RESULT  C  IS  NOT  SET  TO  ZERO 

MATMPY 

c 

MATMPY 

DIMENSION  A(L,J),B(M,K),C(N,K) 

MATMPY 

DO  2  12=1, K 

MATMPY 

DO  2  11  =  1,1 

MATMPY 

C(U,I2)=0. 

MATMPY 

DO  2  13= U 

MATMPY 

2 

C(U,I2)  =  C(U,I2)+A(U,13)*B(I3,I2) 

MATMPY 

RETURN 

MATMPY 

END 

MATMPY 

1 
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